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The dynamics, appearing after a quantum quench, of a trapped, spin-orbit coupled, dilute atomic 
gas is studied. The characteristics of the evolution is greatly influenced by the symmetries of 
the system, and we especially compare evolution for an isotropic Rashba coupling and for an 
anisotropic spin-orbit coupling. As we break the rotational symmetry by making the spin-orbit 
coupling anisotropic, the underlying classical model is chaotic and the quantum dynamics is affected 
accordingly. Within experimentally relevant time-scales and parameters, the system thermalizes in 
a quantum sense. The corresponding equilibration time is found to agree with the Ehrenfest time, 
i.e. we numerically verify a ~ \og(h~ 1 ) scaling. Upon thermalization, we find the equilibrated 
distributions show examples of quantum scars distinguished by accumulation of atomic density for 
certain energies. At shorter time-scales we discuss non-adiabatic effects deriving from the spin-orbit 
coupled induced Dirac point. In the vicinity of the Dirac point, spin fluctuations are large and, even 
at short times, a semi-classical analysis fails. 

PACS numbers: 03.75.Kk, 03.75.Mn 



I. INTRODUCTION 

The physics of ultracold atomic gases has greatly ad- 
vanced in recent years pQ . The high control of system pa- 
rameters, together with the isolation of the system from 
its environment, have made it possible to use such se- 
tups to simulate various theoretical models of condensed 
matter physics [TJ [2] . Of significance in many condensed 
matter models is the response to external magnetic fields. 
Since atoms are neutral, there is no direct way to imple- 
ment a Lorentz force in these systems. Early experiments 
created a synthetic magnetic field via rotation [3] . While 
simple theoretically, these methods are impractical for 
certain setups, and they are limited to weak, uniform 
fields. The first experimental demonstration of laser- 
induced synthetic magnetic fields for neutral atoms [4], 
on the other hand, paves the way for an avenue of new 
situations to be studied in a versatile manner [5HZ] . Ow- 
ing to numerous fundamental applications in the con- 
densed matter community [8] [9], maybe the most im- 
portant direction appears when the laser fields induce 
a synthetic spin-orbit (SO) coupling. Indeed, a certain 
kind of SO-coupling for neutral atoms has already been 
demonstrated [10 , and it is expected that more general 
SO-couplings will be attainable within the very near fu- 
ture [H [12]. 

While SO-couplings can in principle bear identical 
forms in condensed matter and cold atom models, there is 
an inevitable difference, often overlooked, between these 
two systems. The presence of a confining potential for the 
atomic gas can qualitatively change the physics [TJ [3], 
and has only recently been addressed [T3HT7] . Fur- 



* Electronic address: jolarson@fysik.su.se 



thermore, most of these studies are concerned with 
ground/stationary state properties of the system [T3l - 
[T5] , while few works discuss dynamics or non-equilibrium 
physics. Notwithstanding, the experimental isolation of 
these systems suggests that they are well suited for stud- 
ies of closed quantum dynamics [18] . 

Historically, some of the finest experiments regard- 
ing dynamics of closed quantum systems have been per- 
formed in quantum optics p~9j [20]. An early example 
proved quantization of the electromagnetic field by mak- 
ing explicit use of quantum revivals [2T] , Such quantum 
recurrences, in general connected to integrability or small 
system sizes, are now well understood. The situation be- 
comes more complex for non-integrable systems [TS] or 
systems with a large number of degrees-of- freedom [22] . 
One particularly interesting question is whether any ini- 
tial state relaxes to an asymptotic state, and if so, what 
are then the properties of this "equilibrated" state and 
the mechanism behind the equilibration. Both these 
questions have inspired numerous publications during 
the last decade, both theoretical [23j [24] as well as ex- 
perimental [25-27 . A rule of thumb is that in order 
for a closed quantum system to thermalize, i.e. all ex- 
pectation values can be obtained from a microcanonical 
state, its underlying classical Hamiltonian should be non- 
integrable [TS]. While true in most cases studied so far, 
exceptions to this hypothesis has been found [28]. More- 
over, the behavior near the transition from regular to 
chaotic dynamics, classically explained by Kolmogorov- 
Arnold-Moser theory [29 , is not well understood for a 
quantum system [30]. It is therefore desirable to study 
a system where these two regimes can be explored by 
tuning an external parameter, and for which the exper- 
imental methods in terms of preparation and detection 
are already well developed. 

Motivated by the above arguments, in this paper we 
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consider dynamics of a trapped SO-coupled cold dilute 
atomic gas. The SO-coupling is assumed tunable from 
isotropic (Rashba-like) to anisotropic, and hence the sys- 
tem can be tuned between regular and chaotic. Note 
that even though this crossover is generated by a change 
in the form of the SO-coupling, the confining trap causes 
the system to become non-integrable. We distinguish 
between short and long time evolution, where by "long 
time" we mean times similar to the Ehrenfest time. In 
fact, the corresponding time-scale for the thermalization 
is found to agree with the Ehrenfest time, and thereby 
scale as \og(h~ 1 ) / X where A is the maximum Lyaponov 
exponent. This scaling for the thermalization has been 
conjectured in Ref. [3T] . but was not numerically ver- 
ified in these works. At shorter times when the wave 
packet remains localized, we especially study the rapid 
changes in the spin as the wave packet evolves in the 
vicinity of the Dirac point (DP). For energies below the 
DP (£" < 0), we utilize an adiabatic model derived in the 
Born-Oppenheimer approximation (BOA) [32]. Aside 
from some special initial states, we encounter thermal- 
ization in all cases. These exceptions correspond to 
states evolving within a regular "island" in the otherwise 
chaotic sea. Among the thermalized states, the equi- 
librated distributions are found to show quantum scars 
originating from periodic orbits of the underlying classi- 
cal model. The experimental relevance of all our theoret- 
ical predictions are discussed and put in a state-of-the-art 
experimental perspective. 

The paper is outlined as follows. The following sec- 
tion introduces the system Hamiltonian and discusses 
its symmetries. Section [II B| derives the adiabatic model 
by imposing the BOA. A semi-classical analysis, demon- 
strating classical chaos for anisotropic SO-couplings, is 



presented in Sec. Ill The following section considers the 



full quantum model at short times, Sec. IV A and long 



times, Sec. IV B Section IV C contains a discussion re- 
garding experimental relevance of our results. Finally, 
Sec. |V| gives some concluding remarks. 



II. SPIN-ORBIT COUPLED COLD ATOMS 
A. Model spin-orbit Hamiltonian 

Several proposals exist for implementing spin-orbit 
couplings in cold atoms [33U35] . In general, these syn- 
thetic spin-orbit fields are generated through the appli- 
cation of optical and Zeeman fields to produce a set of 
dressed states that are well separated energetically from 
the remaining dressed states [5]. We denote these states 
as pseudo-spin, but emphasize that there is no connection 
to real space rotations. Spatial variation of the dressed 
states will couple the pseudo-spin to the orbital motion 
of the atom. An atom prepared in a pseudo-spin state 
will therefore see an effective Hamiltonian, provided the 
atom is sufficiently cold. 

For a specific configuration of optical fields, one can 



induce the effective Hamiltonian [35] 
p2 i 

H S o = ^ + -muj 2 r 2 + v x p x a x + v y p y cr y , (1) 

where p = (p x ,p y ) is the momentum operator, r = (x,y) 
is the position operator, m is the mass of the atom, and 
lj the frequency of a harmonic trap. The operator &i 
is the i-th Pauli matrix in pseudo-spin space, and the 
velocities V{ couple pseudo-spin to an effective momen- 
tum dependent Zeeman field, B(p) = (v x p x ,v y p y ). This 
momentum-dependent Zeeman field can simulate any 
combination of the Rashba [38] and Dresselhaus [39] SO- 
couplings experienced in semiconductor quantum wells 
and systems alike. 

In the absence of a trap, uj = 0, the spectrum of ([I]) is 



E»(Px,Py) = 7T- (Pi +Pl) + ^xPx) 2 + (VyPy) 2 (2) 
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with the corresponding eigenfunctions 
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is a spinor with helicity \i = ±1 and ip = 
aictaiL(v y p y /v x p x ). These states have well defined mo- 
mentum, but have no velocity since (r) = (V p H) = 0, 
provided the optical fields are maintained. Note further 
that the eigenstates are parametrically dependent on p x 
and p y . 

We remark that for an isotropic SO-coupling, v x =v y , 
the Hamiltonian ([TJ is equivalent to the dual E x e Jahn- 
Teller model, frequently appearing in chemical/molecular 
physics and condensed matter theories [37]. With a 
simple unitary rotation of the Pauli matrices, the SO- 
coupling attains the more familiar Rashba form [38] (or 
equivalently Dresselhaus form [39]). For v x ^ v y , i.e. 
when the SO-coupling is anisotropic, the model becomes 
the dual E x (f3 x + /3 y ) Jahn- Teller model [37 . In par- 
ticular, we have that the total z angular momentum 
J z = L z + is a constant-of-motion for the isotropic 
but not for the anisotropic model. More precisely, break- 
ing of the SO isotropy implies a reduction in symmetry 
from U(l) to Z2. 

Throughout we will use dimensionless parameters 
where the oscillator energy E Q = Tiuj sets the energy- 
scale, I = ^Ti/muj the length-scale, and the characteristic 
time is r = cj _1 . We note that for typical experimental 
setups, uj ~ 10 - 100 Hz and m(v 2 + v 2 )/h ~ 1 - 10 kHz. 
Moreover, in what follows we will refer to pseudo-spin 
simply as spin. When necessary, we introduce a param- 
eter h serving as a dimensionless Planck's constant, i.e. 
hh. In this way, h controls the strength of Planck's con- 
stant and by varying it we can explore how the dynamics 
depends on h. 
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B. Adiabatic model 

The large ratio of the SO energy to trapping energy, 
typically mv 2 /hu - 10 - 1000, suggests that a BOA [32] 
will be valid for experimental implementations. In many 
experiments, the spin will be initially aligned with an 
effective magnetic field. The separation of time-scales 
implies that the spin will adiabatically follow the effective 
magnetic field B(p). We then project into the helicity 
basis to find the adiabatic Hamiltonain, 



H. 



p 2 



ad 



+ V + ? + ? (5) 



V_ + Px 

2 ' 2 2 ' 2 

The trap thus takes the role of kinetic energy and ^ can 
be pictured as a particle in a (dual) adiabatic potential 



fx 



vl 



V,(p x ,p y ) = y + f + ^vlf>l + vffi. (6) 

shown in Fig. [I] for both the isotropic (a) and anisotropic 
(b) cases. In this basis, the "bare states" defined in 
Eq. @ are identified as the adiabatic states [32] . 

We have neglected non-adiabatic corrections arising 
from the vector potential and the Born-Huang term [40] . 
For example, an additional scalar potential 



V nad (PxiPy) 
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will emerge from the action of the SO-coupling on the 
spinor {(f^). This term is order V na( ± ~ (ip\ V 2 \cp) ~ 1/p 2 • 
There will also be an additional vector potential term 
A ~ 1/p. The non-adiabatic corrections diverge near the 
DP, but then fall off rapidly at finite p. The adiabatic 
approximation, i.e. BOA, will be valid if the particle 
avoids p = 0. We will show later that this condition is 
met if the particle is in the lower band and has an energy 
E<0. 

Imposing the BOA, any state propagating on the lower 
adiabatic potential will be denoted &(p x ,p y ,t), and it is 
understood that 



$(Px,Py,t) = </>(p x ,Py,t)\ip-). 



(8) 



The real space wave function fy(x,y,t) is given as 
usual from the Fourier transform of (j)(p x ,p y ,t). 
The time-evolution follows from 4>(Px->Py>t) — 

exp ^— iH^tj (f>(Px,Py: 0). It is also clear that the 

state $>(p x ,p y ,t) determines the spin orientation which 
is inherent in the ket- vector \(fi). More explicitly, the 
time-evolved Bloch vector 



R(t) = (R x (t),R y (t),R z (t)) = ((& x ),(a y ),(a z )) 
takes the form 



(9) 



x(t) = J dp x dp y \(j)(p x ,p y ,t)\ 2 cos(cp), 
(t) = J dp x dp y \(j)(p x ,py,t)\ 2 sm((p), 



y 

Rz(t) 



(10) 
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FIG. 1: Adiabatic potentials of the isotropic (a) and 
anisotropic (b) SO-coupled models. In (a), the lower adi- 
abatic potential V-(p x ,p y ) has the characteristic sombrero 
shape. By considering an anisotropic SO-coupling, the rota- 
tional symmetry is broken and V- (p x ,p y ) possesses two global 
minima at (p x ,p y ) = (0, ±v y ). 



in the BOA, and it is remembered that the parameter cp 
depends on p x and p y . Note that the Bloch vector pre- 
cesses in the equatorial spin x?/-plane. If the wave packet 
$>(p x ,Py,t) is sharply localized, a crude approximation 
for the Bloch vector is given by 



R x (t) = 



Ry{t) 



V x Px(t) 



y '(V^t)) 2 + (VyPy(t)f 
VyPy{t) 



^ {V X P X (t)f + {VyPy(t)f 

R z (t) = o, 

where p a (t) = J dp x dp y \$(p x ,p y ,t)\ 2 p a with a 
III. CLASSICAL DYNAMICS 



(11) 

(12) 

(13) 
x, y. 



Quantum chaos is often defined by having an under- 
lying chaotic classical model. For the full model ([!]), 
the spin degrees-of-freedom cannot be eliminated in a 
straightforward manner in the vicinity of the Dirac point 
and as a consequence it is not a priori clear what the 
underlying classical model would be in this regime. On 
the other hand, in the BOA, the adiabatic Hamiltonian 
can serve as our classical model Hamiltonian. Still, 

it should be noted that we assume (H^) <C 0, such that 
the spectrum contains a sufficiently large number of en- 
ergies below E = 0. Furthermore, we point out that jus- 
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FIG. 2: Poincare sections of the Rashba SO-coupled adiabatic 
model (J5| for the intersections y = (a) and p y = (b). 
The initial energy is E — —192, the SO-coupling strengths 
v x — v y = 30, and the number of simulated semi-classical 
trajectories 18. 



tification of the BOA does not necessarily imply approval 
of a semi-classical approximation which depends on the 
system energy and the actual shape of the dual poten- 
tial V-(p x ,p y ). Nevertheless, as we will demonstrate in 
the following, for the chosen parameters, the agreement 
is indeed very good. 

The classical equations-of-motion of the Hamiltonian 



Vx 
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Vy 



Vx 



Vy 



v lVx 



yjv 2 p 2 +v 2 p 2 

v lVy 
yJ v xVl+VyVl 



(14) 

(15) 
(16) 

(17) 



For the Rashba SO-coupling, v x = v y = v, there is one 
unstable fix point (p x ,p y ) — (0>0) an d a seam of sta- 
ble fix points p 2 + Vy — v<2 , see Fig. [I] (a). For the 
anisotropic case, v y > v X: there are three unstable fix 
points, (p x ,p y ) = (0,0) and (p x ,p y ) = (±v x ,0), while 
there are two stable fix points (p x ,p y ) — (fi^ v y)i see 
Fig. 0(b). 

The classical energy E(x,p x ,y,p y ) = p 2 x /2 + p 2 y /2 + 
x 2 /2 + y 2 1 2 — ^Jv 2 p 2 + v 2 p 2 determines a hypersurface 

in phase space for any given energy E(x,p x ,y,p y ) — Eq. 
The semi-classical trajectories (x{t),p x {t), y(t),p y {t)) live 
on this surface. For the integrable case, v x = v y , 
these surfaces form different tori characteristic for quasi- 
periodic motion. As the rotational symmetry is slightly 
broken, v x ^ v y , the tori deforms and the motion loses 
its quasi-periodic structure [29]. This is the generic 
crossover from regular to chaotic classical dynamics. As 
will be discussed further below, even in the chaotic 
regime, periodic orbits may persist and will greatly af- 
fect the dynamics, both at a classical and a quantum 
level [41]. 

The semi-classical behavior for our system is visualized 
using Poincare sections [42] in Figs. [5] and [3j We solve 
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FIG. 3: Poincare sections of the anisotropic SO-coupled adia- 
batic model |5]) for y = (a) and (c), and for p y = (b) 
and (d). The initial energies are E = —88 (a) and (b), 
and E — —192 (c) and (d), and the SO-coupling strengths 
v x = 20 and v y = 30 for both cases. The corresponding maxi- 
mum Lyaponov exponents have been derived to A « 0.12 and 
A = 0.090 respectively. The number of semi-classical trajec- 
tories is the same as for Fig. [2] namely 18. 



the set of coupled differential equations (14) using the 
Runge-Kutta (4,5) algorithm modified by Gear's method, 
suitable for stiff equations. We have also numerically 
verified our results employing different algorithms [43] . 
In the first figure we display the Poincare sections in the 
xp x plane for the intersections determined by y = (a) or 
p y = (b) of the isotropic model with the SO-coupling 
amplitudes v x = v y = 30. The initial energy is taken 
as E = —192, well below the DP, consistent with the 
BOA. In (b), the section defined by p y = 0, the evolution 
results in ellipses in the Poincare section, characteristic 
of quasi periodic motion. The structure of the Poincare 
section for y = (a) is somewhat more complex. This can 
be understood from the sombrero shape of the adiabatic 
potential V-(p x ,p y ) m , for given x = x' , p x = p' x , y = 0, 
and energy Eq , there are four possible values of p y , and 
this multiplicity of possible p y 's allow the "curves" in 
Fig. [2] (a) to cross. 

Figure [3] presents two examples for anisotropic SO- 
couplings, both with v x = 20 and v y = 30. The 
quasi-periodic evolution is lost and the dynamics become 
mixed, with regions of both chaos and regular dynamics. 
The same conclusions were found in Ref. [44] where a re- 
lated Jahn- Teller model was studied. The two lower plots 
consider the same energies as in Fig. [2] i.e. E = -192, 
while for (a) and (b) E = —88. Expectedly, the higher 
energy increases the accessible volume of phase space. 
For both energies we find islands free from chaotic tra- 
jectories. As will be demonstrated in the next section, 
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within these islands the evolution is regular and the sys- 
tem does not thermalize. The plots also demonstrate 
clear structures also appearing in the chaotic regimes of 
the Poincare sections in which the density of solutions 
changes. 

IV. QUANTUM DYNAMICS 

The idea of this section is to analyze how the corre- 
sponding quantum evolution is affected by whether the 
classical dynamics is regular or chaotic. Of particular im- 
portance is the long time evolution in which the system 
state may or may not equilibrate. However, we study 



also the short time dynamics arising for a localized wave 
packet traversing the Dirac point. In this regime, clearly 
the classical results of the previous section does not hold. 

To study the system beyond the classical approxima- 
tion, we solve the time-dependent Schrodinger equation, 
represented by the Hamiltonians or ([5]), to obtain 
the corresponding wave function ty(x : y,t) at time t. 
Note that for the full model ([!]), the wave function con- 
tains the spin degree-of- freedom ^(x, y, t) = i/)^(x, y,i)\ t 
) + ijj±(x,y,t)\ 1). The non-equilibrium initial state ap- 
pears after a quench in the center of the trap. We prepare 
the system in a quasi-ground state for a shifted trap, and 
at t = suddenly move the trap center to x s = y s = 0, 



V(x,y) 



(x-x s ) 2 {y-y s f 



x s 7^ and/or y s ^ 0, t < 0, 
x s = Vs = 0, t > 0. 



(18) 



By "quasi-ground state" in an anisotropic SO-coupled 
system, we consider an initial state predominantly pop- 
ulated in one of the two minima of the adiabatic po- 
tential V-(p x ,p y ). This seems experimentally reasonable 
where small fluctuations will favor one of the two min- 
ima. For the isotropic case, the phase of <&(p x: p y ,t = 0) 
is taken randomly in agreement with symmetry breaking. 
Given the evolved states ^(x^y, t), we are interested in 
the Bloch vector (10) or its components, and the distri- 
butions \<&(p x ,p y ,tj^ and \^(x,y,t)\ 2 . 

The numerical calculation is performed employing the 
split- operator method [45 which relies on factorizing, for 
short times St, the time-evolution operator into a spatial 
and a momentum part. For small SO-couplings v x and 
v y , the method is relatively fast, while as v x and/or v y are 
increased the time-steps St must be considerably reduced 
and the necessary computational power rises rapidly. In 
addition, for large v x and v y , the grid sizes of position 
and momentum space must be increased, which also in- 
creases the computation time. Thus, we will limit the 
analysis to SO-couplings v x , v y < 30. Furthermore, we 
have found by convergence tests that the full model 
requires much smaller time-steps St than the adiabatic 
one ([5|, and most of our simulations will therefore be re- 
stricted to energies E < for which the BOA is justified. 

The full quantum simulations are complemented by the 
semi-classical truncated Wigner approximation (TWA), 
which has turned out very efficient in order to re- 
produce quantum dynamics [46]. The TWA considers 
a set of N different initial values (xi,yi,p X i,p y i) ran- 
domly drawn from the distributions |^(x, 0) | 2 and 
\&(p x ,Py, 0)| 2 - These are then propagated according to 
the classical equations-of- motion (14). The propagated 



set (xi(t),yi(t),p X i(t),Pyi(t)) gives the semi-classical dis- 
tributions, from which expectation values can be evalu- 



ated. 



A. Short time dynamics 

Before investigating the prospects of thermalization, 
we first consider short time dynamics, by which we mean 
time-scales where the wave packet remains localized. In 
this respect, it is tempting to think of the dynamics as 
semi-classical. However, in the vicinity of the the DP any 
classical description would fail. Equivalently, the spin 
degrees-of-freedom will show large fluctuations which are 
difficult to capture classically. The short time dynamics 
is consequently most interesting for situations with ener- 
gies E > where both the semi-classical approximation 
and the BOA break down, implying that the simulation 
is performed using the full model Hamiltonian 0. For 
these energies, the wave packet can traverse the DP and 
population transfer between the two adiabatic potentials 
V^PxiPy) typically occurs. It is known that such non- 
adiabatic transitions can play important roles for the dy- 
namics, and that the actual transition probabilities be- 
tween the two potentials may be extremely sensitive to 
small fluctuations in the state STJ[48j. In this subsection 
we especially address such non-adiabatic effects. 

There are indeed several relevant time-scales in the dy- 
namics: (i) The spin precession time T sp gives the typical 
time for spin evolution and is proportional to the effective 
magnetic field |B(p)|, (ii) the classical oscillation period 
T c i = 27r, and (Hi) the thermalization time T^, which 
estimates the time it takes for the system to thermalize, 
i.e. when expectation values become approximately time 
independent. Normally, the magnitudes of these times 
follow the list above (in growing order), except in the 
vicinity of the DP where T sp ~ T c \ or even T sp <C T c \ 
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(a) 








u 








FIG. 4: Bloch vector components (dashed lines) and R y 
(solid lines). For the upper plot (a), the trap has been dis- 
placed in th ^/-direction, x s = and y s — 28, while in the lower 
plot (b) the displace direction is the perpendicular, x s = 28 



and y s = 0. In both figures, v x 
average energy E « 280. 
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not spend much time near the DP so the wave packet 
derealization occurs more slowly. To a large extent the 
evolution is driven by harmonicity, in contrast to the ex- 
ample of Fig. [5] (a) where the anharmonicity of the Born- 
Huang term, and the non-adiabatic transitions near the 
DP, push the system away from semi-classical evolution. 
The figure demonstrates how the dynamics can depend 
on the initial conditions, in both (a) and (b), E « 280 but 
the wave packet broadening starts earlier in (a) than in 
(b). This type of state-dependence has been discussed in 
Ref. [47] ; generically there is a period t s where the width 
of the wave packet stays nearly constant, followed by a 
rapid broadening. The time-scale t s depends strongly 
on the initial conditions, while the proceeding evolution 
after t s seems pretty generic for chaotic systems. 



B. Long time dynamics; thermalization 



very close to the DP. While the first two are well defined, 
defining the last one is non-trivial. We can say that (i) 
and (ii) characterizes short time-time scales, and (Hi) 
long time-scales. As will be numerically demonstrated, 
the thermalization time turns out to scale as log(/i _1 )/A, 
where h is the effective dimensionless Planck's constant 
and A the maximum Lyaponov exponent. This suggests 
that the thermalization time agrees with the Ehrenfest 
time 



T E = log(V//i)/A, 



(19) 



with V the effective occupied phase space volume. Te 
is also the typical time-scale where semi-classical (TWA) 
expectation values no longer agree with quantum expec- 
tation values, which can be seen as a breakdown of Ehren- 
fest's theorem [49] . 

From the form of the non-adiabatic coupling 0, it fol- 
lows that transitions between the adiabatic states Q are 
restricted to the vicinity of the DP. These non-adiabatic 
transitions are manifested as rapid changes in the Bloch 
vector ( [To] ). In Fig. [5] we present two examples of the 
Bloch vector evolution (in both examples R z (t) ~ 0). In 
Fig. [5] (a), the trap has been shifted in the ^-direction. 
For short times, the shift of the trap induces a build-up of 
momentum in the opposite ^-direction as a consequence 
of the Ehrenfest theorem. This adds with the non-zero 
^-component of momentum before the quench. The av- 
erage momentum in the x-direction remains zero and as 
a consequence R x (t) ^ 0, see Eq. (11). 



These dynamics change qualitatively if the trap is 
shifted in the x-direction instead of the ^/-direction. For 
sufficiently large shifts of x 8 , the wave packet will set off 
along the adiabatic potentials and encircle the DP. The 
spin dynamics should therefore not display the same type 
of "jumps" that appear when the wave packet traverses 
the DP. Moreover, since the average momentum in the 
x-direction is in general non-zero, R x (t) will also be non- 
zero. The results are demonstrated in Fig. [5](b). Com- 
pared to the first example in (a), the wave packet does 



Whenever we consider an anisotropic SO-coupling, 
v x 7^ v y , from the Figs. [2] and [3] it is clear how the 
adiabatic classical model becomes chaotic. Beyond the 
adiabatic model, it has been shown [50 that the full 
anisotropic model, i.e., E x (f3 x + f3 y ) Jahn- Teller model, 
is chaotic in the sense of level repulsion [51] of eigenen- 
ergies. For the isotropic E x e Jahn- Teller model, on 
the other hand, the level repulsion effect is not as evi- 
dent, however a weak repulsion also in this model signals 
emergence of quantum chaos [52] . 

The goal of this subsection is to study the long time 
dynamics of the system; specifically if equilibration oc- 
curs, and if so, does the equilibrated state mimic a ther- 
mal state. A distinguishing property of thermal states is, 
for example ergodicity, i.e., the distributions \^(x,y,t)\ 2 
and \®(p x ,Py,t)\ 2 spread out over their accessible energy 
shells. Moreover, for a thermally equilibrated state, the 
distributions show seemingly irregular interference struc- 
tures on scales of the order of the Planck cells, which 
normally become even finer in the Wigner quasi distri- 
bution [53-55 . Non-thermalized states, on the contrary, 
typically leave much more regular traces of quantum in- 
terference in their distributions. While such often sym- 
metrical structures are absent for thermalized states, we 
will demonstrate that thermalized distributions may still 
show clear density fluctuations on scales larger than the 
Planck cells. These are examples of quantum scars and 
they are remnants of classical periodic orbits [4T] . 

We begin by considering the adiabatic isotropic model 
with v x = v y = 30, and trap shifts x s = y s = 16. Af- 
ter a quench of the trap position, the initial energy is 
E = (H^) « —192. This energy corresponds to the 
energy of the Poincare section presented in Fig. [2] The 
resulting distributions are shown in Fig. [6] (a) and (b) 
after a propagation time tf = 400 . The final time t f ap- 
proximates 60 classical oscillations. Both the real space 
density \Sfr(x, y, t)\ 2 and momentum density \<&(p X: p y: t)\ 2 
reveal clear interference patterns as anticipated. The DP 
at the origin (p x ,p y ) = (0,0) repels the wave function 
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FIG. 5: (Color online) Distributions |^(x, £/)| 2 ((a) and 
(c)) and \<$>(p x ,p y ,tf)\ 2 ((b) and (d)) at t f = 400 for the 
Rashba SO-coupled model. At time t — 0, the trap is sud- 
denly displaced from xo — yo — 16 to xo — yo =0. The 
initial ground state is then quenched into a localized excited 
state. The upper two plots (a) and (b) display the results 
from full quantum simulations of the adiabatic model (J5|, 
while the lower plots (c) and (d) show the corresponding 
semi-classical TWA distributions. The average semi-classical 
energy E « —192 with a standard deviation 5E « 22. The 
dimensionless SO-coupling strengths v x — v y — 30. 
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FIG. 6: (Color online) Same as Fig.|5]but for the anisotropic 
SO-coupled model with v x = 20 and v y = 30. The largely 
populated regions are so called quantum scars and derive from 
properties of the underlying classical model, i.e. they are not 
outcomes of some coherent quantum mechanism. 



FIG. 7: (Color online) Same as Fig.|6]but for an initial energy 
E > 0. The dimensionless SO-couplings v x = 14 and v y = 21, 
while the shifts x s = y s = 16 giving an average energy E = 
(# so )~36.5. 



forming a "hole." The lack of zero momentum states in- 
duces a mass flow in real space and a similar "hole" in 
its distribution. The classically energetically accessible 
regions are given by 



x 2 + y 2 < 2E n 



(20) 



pi+p 2 y -2Jv*p>+vffi<2E u 



where E max is the maximum energy component notice- 
ably populated by the state. 

The quantum results are compared with the TWA dis- 
tributions displayed in the lower plots (c) and (d) of the 
same Fig. [6] The same kind of ring-shape is obtained, 
and the concentration in density appears at the same lo- 
cations for both the quantum and classical simulations. 
Expectedly, the quantum interference taking place within 
the wave packet is not captured by the TWA. This fol- 
lows since single semi-classical trajectories are treated 
independently, i.e. added incoherently, while a quantum 
wave packet must be considered as one entity. For a TWA 
approach of the full isotropic E x e Jahn- Teller model ([!]) 
we refer to Ref. [56] . 

The situation is drastically changed when we break 
the rotational U(l) symmetry by assuming v x ^ v y . The 
result for low initial energy is depicted in Fig. [6] (a) and 
(b). The energy is comparable to the potential barrier 
separating the two minima in the adiabatic potentials, 
and as a consequence, the wave packet is predominantly 
localized in the left minima. The density modulations 
seems now much more irregular in comparison to Fig. [6] 
In the seemingly random density distribution, some clear 
density maxima emerge, both in momentum as well as 
in real space. These density accumulations derive from 
periodic orbitals of the underlying classical model and are 
termed quantum scars [4TJ [57] . The appearance of scars 
is an example of the classically chaotic model leaving a 
trace in its quantum counterpart. The scars are also 
captured in the semi-classical TWA, shown in Fig. [6] (c) 
and (d), supporting their classical origin. 

When we shift the trap for larger values on x s and 
y s , the energy is increased and at some point the BOA 
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FIG. 8: Examples of the phase space area A x (t) for different 
h- values (h — 1, 2, 3, . . . , 10). The upper plot (a) gives A x (t) 
without shifting the time, while for the lower one (b) time has 
been shifted by 5 = log(/i)/A. The arrow indicates increasing 
h- values. It is clear how the spread in A x (t) between different 
h values is suppressed when we shift the time. The trap shifts 
x s — y.s = 19 resulting in an energy E w —88. The maximum 
Lyaponov exponent A = 0.18. 



breaks down. An example, obtained from integrating the 
full model 0, is presented in Fig. [7] For these higher 
energies there are no signs of quantum scars. As for the 
situation of Fig. [6j the spread of the wave packet and the 
irregular interference patterns indicates thermalization. 

This far we have demonstrated thermalization for the 
anisotropic SO coupled model, but not discussed corre- 
sponding time-scales. One related question is how the 
evolution of various expectation values scale with h (di- 
mensionless Planck constant). It has been argued that 
the Ehrenfest time, Eq. (19), can be a measure of the 
thermalization time [3Tj . We will now explore how the 
phase space area A a (t) = AaAp a (a = x,y), where 



Aa and Ap a are the variances of a and p a respectively, 
evolves for different values of h. Since A x (t) and A y (t) 
behave similarly we focus only on A x (t). For thermaliza- 
tion, A x (t)A y (t) is an effective measure of the covered 
phase space volume, and for large times t it should more 
or less approach the accessible phase space volume as the 
distribution spreads over the whole energy shell. We have 
chosen to study A x (t) since it fluctuates relatively little 
before reaching its asymptotic value. In Fig. [8] (a) we dis- 
play A x (t) for 10 different values on h ranging from h = 1 
to h = 10. The arrow in the plot shows the direction of 
increasing ft's. As is seen, by increasing h the wave packet 
broadening starts earlier and the state equilibrates faster. 
If the Ehrenfest time T# sets the typical time scale in the 
process, by shifting the time with 5 = \og(h)\ we should 
recover a "clustering" of the curves. This is indeed ver- 
ified in Fig. [8] (b) where the curves have been shifted in 
time by S. The corresponding Lyaponov exponent A has 
been optimized in order to minimize the spread in the 
curves. The obtained value A = 0.18 is somewhat larger 
than the numerically calculated one A = 0.12 but still of 
the same order. The picture also makes clear that the 
wave packet broadening kicks in after some time t s as 




FIG. 9: (Color online) Sections of \ip(x,y = 0)| for different 
values on the dimensionless Planck's constant h: h = 1 (black 
solid line), h = 2 (blue dotted line), and h — 3 (red dashed 
line). The final time tf = 80, x s — y s — 16, and v x = 14 
and v y = 20. As a comparison between classical and quan- 
tum results, we also include the TWA results as a green solid 
line, calculated for h = 1. The green line has been shifted 
downward with 0.02 for clarity. 



anticipated above. 

The route to thermalization can typically be di- 
vided into; (i) a classical drift, and (ii) quantum dif- 
fusion [3Tj . The role of the quantum diffusion for ther- 
malization was analyzed in Ref. [31], where it was found 
to "smoothen" the phase space distributions preventing 
sub-Planck structures. For the classical drift there is no 
lower bound on the fineness of density structures that can 
form, and characteristic for classical chaotic dynamics is 
that ever finer formations build-up as a result of the typ- 
ical "stretching-and-folding" mechanism. However, in a 
quantum chaotic system, when the structures reach the 
Planck cell regime, the quantum pressure becomes too 
strong and the quantum diffusion then prevents any fur- 
ther structures to form. Thus, Planck's constant sets 
a lower bound for the fluctuations in the distributions. 
This quantum smoothening is demonstrated in Fig. [9j 
where we plot a section of \ip(x,y = 0)| for different 
values on the scaled dimensionless Planck's constant h 
(= 1, 2, 3 for black, blue, and red lines respectively). 
The effect is clearly seen in the figure. A similar pattern 
is found (not shown) also for the momentum distribu- 
tions. For the classical system, corresponding to h = 0, 
there is no lower limit on how fine the structures can be. 
We indicate this by also plotting the TWA results in the 
same figure as a green line (note that the green line has 
been shifted downward in order to separate it from the 
quantum results). The number of trajectories used for 
the figure is 250 000, and if we would like to produce 
finer structures (by propagating the system for longer 
times) we would need many more trajectories and the 
simulation would rapidly become very time consuming. 

Related to the above discussion a note on quantum 
phase space distributions is in order. It is well known 
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FIG. 10: (Color online) Same as Fig. [6] but for the shifts 
x s = 20 and y s = 0. For the given dimensionless parameters, 
the initial state is such that its dynamics should be regular 
according to the corresponding Poincare section, Fig. [3] The 
energy E « —250. 



that sub-Planck structures are common in the Wigner 
distribution [53]. This is not contradicting any quan- 
tum uncertainty relation. After all, the Wigner distribu- 
tion is not a proper probability distribution, despite the 
fact that its marginal distributions reproduce the cor- 
rect real and momentum space probability distributions. 
The Husimi Q-function, while not possessing the proper 
marginal distributions, is positive definite and lacking 
singularities, and it is indeed found that the Q-function 
does not support sub-Planck structures [59] . 

We finish this subsection by analyzing the dynamics in 
the islands of the Poincare sections of Fig. [3] where the 
classical theory predicts regular evolution. From Fig. [3] 
(c) we have that for p x « 20 and x « y w the evolu- 
tion should be regular. We can achieve such a situation 
by using the quench-shifts x s — 20 and y s = 0. As for 
the examples above, we propagate the state for a time 
tf = 400, and the resulting distributions are given in 
Fig. [To] The striking difference with Figs. [6] and [7] is evi- 
dent; no irregular structure is apparent, but clear regular 
interference patterns are. We have verified that the in- 
terference structure prevails also after doubling the time, 
t f = 800. 



C. Proposed experimental realization 

Much of the above dynamics can be observed in a sys- 
tem of cold atoms with synthetic SO-coupling, for exam- 
ple, a system of 87 Rb with a synthetic field induced by 
the 4- level scheme [33]. In this system, the recoil energy 
E r = mv 2 ~ h x 50 kHz. The synthetic field limits the 
lifetime of the experiment to t\ ~ Is [3J [10]. To push 
the experiment into the long time regime, we will use 
a trapping frequency of uj/2tt = 30 Hz. These parame- 
ters will give a dimensionless value of v y — ~ 11, 
with v x tunable between and 11. The large trapping 
frequency will provide a sufficient number of oscillations 
for thermalization to occur. We could consider values of 
v y ~ 30 by decreasing the trapping frequency to 10Hz, 
but then the lifetime of the system may be at the boarder 



for thermalization. 

The condensate can be adiabatically loaded to one of 
the two states at the bottom of the momentum-space 
potential, defined by p = ±mv y y. The quench can then 
be preformed by shifting the minimum of the real-space 
trapping potential. We then let the system evolve until 
we reach either the thermalization time, or the lifetime 
of the experiment. The momentum distribution can be 
measured with a destructive time-of-flight (TOF) mea- 
surement [4] [10] , which should reveal thermalization as 
well as signatures of quantum scars. Repeated experi- 
mental measurements allow for time-resolved calculation 
of expectation values. Similarly, the quantum spin jumps 
near the DP, as discussed in Sec. |IV A[ can be observed 
using a spin-resolved TOF measurement. 

As a final remark, for a weakly interacting gas we 
work near a Feschbach resonance [60]. However, for re- 
alistic parameters [61], we estimate a scattering length 
a s ~ 3 x 10 -9 m, N ~ 5 x 10 5 atoms, and a transverse 
harmonic trapping frequency uj z ~ 100 Hz. For these 
parameters, the characteristic scale of the non-linearity 
is /jl ~ h x IkHz, which is smaller than the recoil energy 
above, suggesting the non-linear term will play only a 
minor role. We have numerically verified that the results 
do not change qualitatively by solving the corresponding 
non-linear Gross-Pitaevskii equation. Indeed, we find the 
deviations with a non-linearity are not large enough to 
be seen by eye. 



V. CONCLUSIONS 

In this paper we studied dynamics, deriving from a 
quantum quench, in anisotropic SO-coupled cold gases, 
focusing primary on aspects arising from the fact that 
the underlying classical model is chaotic. The evolution 
of the initially localized wave packet on its way to equili- 
bration has been analyzed, and we have shown how a clas- 
sical period of limited spreading is followed by a collapse 
regime dominated by rapid spreading. After the collapse 
period, the wave packet is maximally delocalized, but 
still possesses quantum interference structures. At the 
Ehrenfest time, the state has approximately equilibrated 
as is seen in the decay of expectation values, as well as 
seemingly irregular density fluctuations both in real and 
momentum space. We showed that the fine structure of 
these fluctuations are limited by the quantum diffusion, 
and thereby the size of the Planck's constant h. For the 
isotropic model, after the collapse no thermalization is 
found, as is expected from the integrability of the under- 
lying classical model. 

For smaller energies, when the wave packet predomi- 
nantly populates one of the dual potential wells, thermal- 
ization is again seen. Here, however, an additional phe- 
nomenon appears in terms of quantum scars. These den- 
sity enhancements emerge along classically periodic or- 
bits. They are classical in nature and long lived. We also 
demonstrated that for certain fine tuned initial states, 
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the dynamics stays regular even in the anisotropic model. 
In the classical picture, these solutions correspond to the 
ones belonging to regular islands in the otherwise chaotic 
Poincare sections. 

We argue that the present system is ideal for studies 
of quantum chaos and quantum themralization for nu- 
merous reasons. The system parameters can be tuned 
externally by adjusting the wavelength of the lasers in- 
ducing the SO-coupling, and as we discussed in Sec.|IVC 



the SO dominated regime is reachable in current exper- 
iments. Moreover, both state preparation and detection 
are relatively easily performed in these setups. Equally 
important, the system is well isolated from any environ- 
ment and coherent dynamics can be established up to 
hundreds of oscillations which is well beyond the themral- 
ization time. The energy of the state is simply controlled 
by the trap displacement, and it should for example be 
possible to give the system small energies such that the 
atoms reside mainly in one potential well where quantum 
scars develop. 

We finish by pointing out that the present model is also 
different from most earlier studies on quantum thermal- 



ization p~8l |24] in the sense that the dynamics is essen- 
tially "single-particle" and not arising from many-body 
physics. Related to this, we have numerically verified 
that adding a non-linear term g\^(x,y,t)\ 2 to the Hamil- 
tonian does not change our results qualitatively for mod- 
erate interaction strengths g. 
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